% KurmannOtrok_main.m

% Main program to replicate results in
%
%   Kurmann, A. and C. Otrok (2012). "News Shocks and the Slope of the Term Structure of Interest rates." American Economic Review (forthcoming)
%
% See readme.txt file for an explanation of the different programs.
%
% Written by Andre Kurmann, Federal Reserve Board. Last version: December 2012
% 
% Some small utilities used in the code are from J.P. Lesage's econometrics
% toolbox and are clearly acknowledged.
%
% Use of code for research purposes is permitted as long as proper reference to source is given. 
%
%------------------------------------------------------------------------------------------------------------------------------------------------

%clear all
%close all
%clear p nvars nlags p data 

%-------------------------------------------------------------
% Step 0: parameter settings
%-------------------------------------------------------------
nlags = 4;                    % number of lags in VAR ; originally =4  
nconst = 0;                   % 1 = constant for , 0 otherwise :Even with no constant term estimation, 
                              % ** CHOOSE nconst=1 since we shut down a constant term in "three_var_final_hiro", "olsvarc_hiro" and related files 
                              % ** by forcing the estimate of a constant term to be zero (V matrics to be zero in "olsvarc_hiro").  
                             
prior=1;                      % for prior,  1 =  standard prior
                              %             inf = flat prior
ndraws=0;                  % number of draws when computing error bands
bound=0.165;                   % lower bound of confidence interval
                              % (e.g. for bound=0.16, CI is 16% - 84%)

nimp = 20;                    % horizon of VD and IRFs to be shown                        

KLbar=0;                      % lower horizon bound for max FEV share objective
KUbar=40;                     % higher horizon bound for max FEV share objective        
                                       
est_method = 0; % 0 : OLS point estimation
                % 1 : Bayesian estimation with Minnesota prior 
                % 2 : OLS point estimation and classical bootstrap             
                
shock_extract = 2; % 1: Uhlig's max FEV share shock
                   % 2: TFP news shock
                   % 3: Current TFP shock
                   
slope = 0;  % 1: maximize FEV share of slope if Uhlig method is used
            % 0: maximize FEV share of level if Uhlig method is used
              
use_yields = 0; % 0: use no yield curve variables
                % 1: use spread between long and 3-month rate and long-rate
                % 2: use yield spread b/w long and ffr -> long rate is not
                % in VAR but computed from spread and ffr thereafter
                % 3: use Diebold-Li slope and level factor

%-----------------------------------------------------------------------
% Step 1: load and transform data
%-----------------------------------------------------------------------
esty1=1957;     % First Year for estimation     (max sample: 1959:2 - 2005:2 except if 10-year yield is used in which case sample starts in 1971:4)
estq1=2;        % First quarter for estimation  
%esty2=1999;     % Last Year for estimation      
esty2=2004;     % Last Year for estimation      
estq2=4;        % Last quarter for estimation 
%three_var_final_hiro_var;
%three_var_final_hiro_var;
[yk, xk, varsk, levk]=dataset_quarterly_sw_var(esty1,estq1,esty2,estq2,nlags,shock_extract,use_yields,slope);
%y is T x nvars matrix of rhs variables
%x is T x nvars*nlags of lhs variables
[T,nvars]=size(yk);



%---------------------------------------------------------------------
% Step 3: Estimate VAR, extract shocks, FEV shares explained, and IRFs 
%---------------------------------------------------------------------
if nconst==1;
   const = ones(T,1); 
   xk = [xk const];
   ncoeffs = nvars*nlags+1;
else
   ncoeffs = nvars*nlags; 
end;

%ncoeffs + 1 x nvars matrix of regression coefficients
bb=xk\yk; %OLS coefficients: ncoeffs x nvars (** more efficient way of computing least squares; i.e. b=inv(x'*x)*x'*y;

%three_var_final_hiro_var;

%b=beta_ols;
b=bb;
res = yk - xk*b;          %T x nvar matrix of residuals
vmat = (1/T)*(res'*res);
stds=sqrt(diag(vmat));

%OLS point estimates
if est_method==0; 
    if (shock_extract == 1)
        [data,v]=uhlig(b,res,vmat,nvars,nlags,KLbar,KUbar,nimp,slope,use_yields);
    elseif (shock_extract == 2)
        [data_k,v]=tfpnews(b,res,vmat,nvars,nlags,KLbar,KUbar,nimp,slope,use_yields);
    elseif (shock_extract == 3)
        [data,v]=current_tfp(b,res,vmat,nvars,nlags,nimp,slope,use_yields);
    end
    
%Bayesian estimation with Minnesota prior
elseif est_method==1;
     
    [hm,bm] = minneprc(yk,xk,nlags,1,nconst,levk,prior); %Minnesota prior
    xxx = inv(kron(eye(nvars),xk'*xk) + diagrv(eye(size(hm,1)),hm));
    bb = xxx * (vec(xk'*yk) + bm);    %stacked posterior means of regression coeffs: vec(bpost) = inv[inv(H) + kron(I,X'X)] * [vec[X'Y] + inv(H)*bprior]
    b=reshape(bb,ncoeffs,nvars);    %posterior coefficient matrix 
    res = yk - xk*b;                  %T x nvar matrix of residuals
    vmat = (1/T)*(res'*res);
    
    if ndraws==0;
        if (shock_extract == 1)
            [data,v]=uhlig(b,res,vmat,nvars,nlags,KLbar,KUbar,nimp,slope,use_yields);
        elseif (shock_extract == 2)    
            [data_k,v]=tfpnews(b,res,vmat,nvars,nlags,KLbar,KUbar,nimp,slope,use_yields);
        elseif (shock_extract == 3)
            [data,v]=current_tfp(b,res,vmat,nvars,nlags,nimp,slope,use_yields);
        end

    else
        sxx=chol(xxx)';
        sinv=chol(inv(vmat));
        datadraws=[];
        vdraws=[];
        randn('state',0);
        randmatrix=randn(nvars,T,ndraws);
        randvec=randn(nvars*ncoeffs,ndraws);
        for jjj=1:ndraws;
            [bbdraw,vmatdraw] = niw(b,sxx,sinv,T,randmatrix(:,:,jjj),randvec(:,jjj)); %drawing from posterior coefficient matrix
            bdraw=reshape(bbdraw,ncoeffs,nvars);
            res = yk - xk*bdraw;
            %compute results for each draw
            if (shock_extract == 1)
                [output,v]=uhlig(bdraw,res,vmatdraw,nvars,nlags,KLbar,KUbar,nimp,slope,use_yields);
            elseif (shock_extract == 2)              
                [output,v]=tfpnews(bdraw,res,vmatdraw,nvars,nlags,KLbar,KUbar,nimp,slope,use_yields);
            elseif (shock_extract == 3)
                [output,v]=current_tfp(bdraw,res,vmatdraw,nvars,nlags,nimp,slope,use_yields);
            end 
            datadraws(:,:,jjj) = output;
            vdraws(:,jjj) = v;
        end
        %compute median and confidence intervals   
        datadraws=sort(datadraws,3);
        low=round(ndraws*bound); high=round(ndraws*(1-bound));
        datalow=datadraws(:,:,low);
        datamed=median(datadraws,3);
        %datamed=mean(datadraws,3);
        datahigh=datadraws(:,:,high);
        data(:,:,1)=datamed; data(:,:,2)=datalow; data(:,:,3)=datahigh;
        vmed = median(vdraws,2);   %median values of shock
    end

%OLS and classical bootstrap for CIs       
elseif est_method==2;
    datadraws=[];
    naccept=0;    %counter for number of accepted sdraws
    for j=1:ndraws;
       % bnew=bootb(b,x,res,nlags,nconst);
        bnew=beta_boots(:,:,j);
        
        %bnew_=bootb_hiro(b,x,res,nlags,nconst);%hiro
        %BIAS_=[BIAS V]';%hiro
        %BIAS_=BIAS_(:,1:q);%hiro
        %bnew=bnew_+ BIAS_;%hiro
        
        res = yk - xk*bnew;          
        vmat = (1/T)*(res'*res);
        %compute results for each draw       
        if (shock_extract == 1)
            [output,v]=uhlig(bnew,res,vmat,nvars,nlags,KLbar,KUbar,nimp,slope,use_yields);
        elseif (shock_extract == 2)
            [output,v]=tfpnews(bnew,res,vmat,nvars,nlags,KLbar,KUbar,nimp,slope,use_yields);
        elseif (shock_extract == 3)
            [output,v]=current_tfp(bnew,res,vmat,nvars,nlags,nimp,slope,use_yields);
        end
            datadraws(:,:,j) = output;
            vdraws(:,j) = v;
        end
        %compute median and confidence intervals   
        datadraws=sort(datadraws,3);
        low=round(ndraws*bound); high=round(ndraws*(1-bound));
        datalow=datadraws(:,:,low);
        datamed=median(datadraws,3);
        %datamed=mean(datadraws,3);
        datahigh=datadraws(:,:,high);
        data(:,:,1)=datamed; data(:,:,2)=datalow; data(:,:,3)=datahigh;
        vmed = median(vdraws,2);   %median values of shock
end

%----------------------------------------------------------------------
% Step 3: Plots
%----------------------------------------------------------------------

%FEV shares explained by shock
%plotFEV(data,vars,slope,use_yields,shock_extract);
   
%IRFs to shock
%plotIRF(data,vars,slope,use_yields,shock_extract)

% decomposition of yield curve response into EH part and term premia part
if use_yields == 1 | use_yields == 2
    plotEH(data,varsk,slope,use_yields,shock_extract)
end